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High order second derivative methods 
with Runge—Kutta stability for the 
numerical solution of stiff ODEs 
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Abstract 


We describe the construction of second derivative general linear methods 
(SGLMs) of orders five and six. We will aim for methods which are A-stable 
and have Runge-Kutta stability property. Some numerical results are given 
to show the efficiency of the constructed methods in solving stiff initial value 
problems. 
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1 Introduction 


In many fields such as control theory, chemical kinetics, biology and the 
movement of stars in galaxies, dynamic behavior is modeled by systems of 
ordinary differential equations (ODEs). We consider the autonomous ODEs 
in the form 


y(t) = f(y(@)), © € [20,7], 


1 
y(Xo) = Yo, oe 


where f : R™ — R™ and m is the dimensionality of the system. We restrict 
our attention to autonomous systems because non-autonomous systems can 
be made autonomous by adding an extra equation to the system. 

For system (1), let g := f,f. For problems in which g can be calcu- 
lated along with f, at a moderate additional cost, second derivative meth- 
ods become feasible. General linear methods (GLMs) [6,7, 12] as a unifying 
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framework for the traditional methods, like Runge-Kutta methods, linear 
multistep methods, predictor—corrector methods and hybrid methods, have 
been extended to the second derivative general linear methods (SGLMs) by 
Butcher and Hojjati [8]. These methods which are s-stage and r-value, for 
the numerical solution of (1) are given by 


Yi"! = h(A@ In) FV) + 2?(4@ Im)G(Y") + U @ Im)yl4), 


= (2) 
yl" = h(B @In)F(Y™) + h?(B @In)G(Y!)) + (V @In)y"-4, 


where h is the stepsize, A, A € R°**, U € R°*", B, B € R"™™* and V € 
R’*” and notation @ is the Kronecker product. Here, Y!"! = (yen ys_, is an 
approximation of stage order q to the vector y(@p,_1+ch) = [y(@n_-1+¢;h)]#_1, 
Le. 


Ge ok 
vil = S>S pky® (ey) +O(R), 151,200.58, (3) 
k=0 ~ 
FVM) := [f"))x1 and G(V")) -= [g(¥f"))]=_1 where the vector ¢ = 
[c1 C2 +++ cs] is the abscissa vector. Also the vectors yl"—! = yer yr, 
and yl"] = yr, are the input and output vectors at the step number n, 
respectively, which for a method of order p take the following forms 
P 
yee _ os igh y™ (fn—1) + One?) i= 1,2,...,7, (4) 
k=0 
and . 
B= Vrach y (en) + OW), FH 5) 
k=0 


for some a; € R associated with the method. 

The main features of SGLMs including pre-consistency, consistency, zero- 
stability and types of these methods have been discussed in [3]. It has been 
shown in [4] that the SGLM (2) with the input vector (4) has order p and 
stage order q = p iff 


e* = zAe® + z7Ae* + Uw(z) + O(2?*7), (6) 
e*w(z) = zBe® + 27 Be® + Vu(z) + O(2?*1). (7) 
where 
se 
ee = lew 02? ele? 


and w(z) is a vector with elements given by 


Pp 
k 2 
wi(z) = ) QiRe”, t=1,2,---,r. 
k=0 


High order second derivative methods with Runge-Kutta ... 3 


In the special SGLMs with p q r 8s, U I, and Ve =e,e = 
[1,1,...,1]" © R*, an equivalent condition for order conditions has been 
found in [5] as 


B= Bo— AB, — AB2 — V B3 — (B—VA)B4+ VA, 
where the (i, 7) elements of Bo, Bi, Bo, B3, and By are given respectively by 


Io" oi(a)de—gi(l+ei) (1+) 0 Gi(a)de ——-$5(i) 
OG). ° Oe) ° bj (ej) ” bj(e;)” bj (cj) 


Here, 


s 
gi(x) = II (x — c;), i=1,2,---,s. 
J=1,j78 

Construction of SGLMs which are also suitable for the numerical solution of 
differential algebraic equations (DAEs) has been discussed in [10]. Some ob- 
tained order barries for different types of SGLMs, found in [3,4,10], are useful 
in construction of these methods. These barriers have been also confirmed 
by means of order arrows by Abdi and Butcher [1,2]. Recently, efficiency of 
these methods in solving stiff ODEs arising from chemical reactions has been 
shown in [11]. 

In continuation of studying on SGLMs, in this paper we construct A- 
stable methods of orders five and six with r = s = 3 and Runge-Kutta 
stability property. 

Next sections of this paper are organized as follows: In Sec. 2, we discuss 
about stability behaviour of Runge-Kutta stable three-stage methods. Sec. 
3 is devoted to construction of SGLMs of orders five and six with A-stability 
property. Some numerical experiments are given in Sec. 4 to demonstrate 
the efficiency of the constructed methods. 


2 RKS three-stage methods 


We first recall that the stability matrix for SGLMs can be obtained by ap- 
plying the methods to the standard test problem of Dahlquist [9] y’ = Cy, 
where ¢ is a complex number, which it is 


M(z)=V + (2B+ 2B)(I-zA- 2A) 'U, 


where z = h¢. Thus, we are interested in stable behavior of powers of M(z). 
If M(z) has only a single non-zero eigenvalue, R(z), then the method is said 
to possess Runge-Kutta stability (RKS) property. For RKS methods, the 
stability behaviour is related to R(z). 

For the methods in which coefficient matrices A and A are lower triangular 
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with the same elements A and yw on the diagonal, respectively, R(z) takes the 
form 


R(z) = wooo (8) 


where deg(N) < 2s. For the methods of order five and six with three stages 
that will be discussed in Sec. 3, the polynomial N defined in (8) satisfies 


N(z) = (1— Az — p2?)8e? — C5z® + O(2"), 
and 
N(z) = (1— Az — pz?)3e7 + O(2"), 
respectively, for an arbitrary C's as the error constant of the method. For the 


method of order six, the error constant is 


1 1 1 


_ : 1 3 
~ 5040 240 40 


2) 22 


1 1 1 
C6 Le put ( ah)» + ( 


40 24 


For these methods to be A-stable, using E-polynomial theorem [7], it is 
necessary and sufficient that A > 0, 4 < 0, and so that the E(y) is non- 
negative for y real where the E-polynomial is defined by 


E(y) =| — Aiy + py?|° — |N(iy)?, 
where i is the imaginary unit. The boundary of the regions of A-stable choices 


of (A, 4) for the methods of order five (with different values of C5) and order 
six are plotted in Figure 1 and Figure 2. 


0.2 


0.4 


—0.6 


0 0.5 1 1.5 2 


m 
Figure 1: The boundary of the regions of A-stable choices of (A, ) for s = 3, p = 5 
corresponding to C = —10~3, —5 x 10-4, -10~4 
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Figure 2: The boundary of the region of A-stable choices of (A, 4) for s = 3, p=6 


3 A-stable RKS methods of orders 5 and 6 


Construction SGLMs of orders p = q < 4 has been discussed in [3-5, 10]. In 
this section, we construct A-stable three-stage methods of orders five and 
six with RKS property. Throughout the construction of these methods, we 
will consider U = I, and V = ev? where v € R” and v¥e = 1. The later 
guarantees zero-stability of the methods [3]. 


3.1 Order 5 methods 


Choosing c = [0 $ 1J”, (Au) = (0.6,—0.1) from the intersection of the 


regions in Figurel and solving the order conditions and the nonlinear RKS 
conditions, the coefficients matrices of the method take the following forms 
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0.6000000000 0 0 
0.4538633794 0.6000000000 0 
0.8442059328 0.8999163314 0.6000000000 


—0.1000000000 0 0 
—0.1450566118 —0.1000000000 0 
—0.9847293116 —0.1278647721 —0.1000000000 


A= 


0.3902646263 0.4639576064 0.2524239604 
—0.3312778090  1.1306242731 0.3534363496 
5.0478598121 —4.1644469839 —0.5208888994 


—0.2677332867 —0.3732899225 —0.0223237563 
—0.4095181371 —0.6362626571 —0.0357186615 
0.5750983052  1.6053219094 0.0622616286 


B= 


[1. 2203054517 —0.3423946125 0. 1220891608] ” . 


This method is A-stable with the error constant Cs % —3.50 x 107¢. 


3.2 Order 6 methods 


Choosing c = [0 c; 1]”, c, as a free parameter, and solving the order condi- 


tions and the nonlinear RKS conditions, we get c; = —1.4989329045 and the 
coefficients matrices of the method take the following forms 


[° 4007120047 0 0 
A= |‘ 5974459850 0.4007120047 0 ; 
0.7281456081 0.0121320319 0.4007120047 
—0.0612701047 0 0 
A= | —0.0145743957 —0.0612701047 0 ; 
0.3881180321 0.1117302066 —0.0612701047 
| 1.1371686053 0.2249968367 0.0903218055 
B= | —0.0512895056 0.1078326109 —0.6604347472 | , 
1.5642870990 0.3929237249 —0.2450012162 
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—0.0425486219 0.0078897842 —0.0128566928 
BS 0.1945434509 —0.0296649869 0.0449770864 
0.3584398092 0.0701030286 —0.0116769898 


y] 


v= [ 0.8572479903 0.2113738061 —0.0686217964]” . 


The obtained value for (A, js) is the interior of the region of A-stable choices 
presented in Figure 2. The error constant for this A-stable method is Cg = 
2.56°« 107°. 


4 Numerical verifications 


In this section we present some numerical results by applying the constructed 
methods of orders five and six in Sec. 3, in order to demonstrate the theo- 
retical expectations. Computational experiments are carried out by applying 
the methods to the following two stiff problems. 


S1— The non-linear stiff test problem 
yi (x) = —1002y; (x) + 1000y3 (a), yi(0) = 1, 
y9(x) = y1(x) — yo(x)(1 + yo(x)), yo(0) =1 


The exact solution is y;(#) = exp(—2x) and yo(%) = exp(—x) and 
x € (0, 1]. 


S2— The stiff initial value problem arose from a chemistry problem 


y} (x) = —0.013y2 — 1000y1y2 — 2500y1y3, yi (0) = 0, 
y(x) = —0.013y2 — 1000y1 2, y2(0) = 1, 
y3(x) = —2500y1 ys, y3(0) = 1. 


The reference solution at x = 2 is 


yi(2) = —0.3616933169289 x 10~°, 
y2(2) = 0.9815029948230, 
y3(2) = 1.018493388244. 


Numerical results for the Problem $1, reported in Table 1, illustrate accuracy 
of the methods of order 5 and 6. These results are obtained with fixed 
stepsizes h = 1/2* with several integer values for k. In this table, we have 
listed norm of error |le;(x)|| at the endpoint of integration « = 1. Also, 
in this table, the rows p refer to the numerical estimates to the order of 
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convergence, computed by the formula p = logs(||en(«)||/||en/2(x)||) where 
en (x) and ep 2(x) are errors corresponding to stepsizes h and h/2. 


Table 1: The global error at the end of the interval of integration [0,1] for 
problem S1 

h O72 2-3 2-4 2-5 

Order 5 method | 2.25x 107" 5.611079 1.51x1071!9 4.34 x 10-12 


p 5.33 5.22 5.12 


Order 6 method | 6.92 x 1078 2.94x107!9 2.45x107!2 5.03x107!4 


Pp 7.88 6.91 5.61 


Numerical results for the Problem $2 are given in Table 2 with stepsize 
h = 0.001. Comparing the obtained results by the methods with the refer- 
ence solution shows the efficiency of the methods for solving stiff non-linear 
problems. 


Table 2: Numerical results for problem $2 solved by the methods of orders 
five and six 


x y Order 5 method Order 6 method 

YL —0.3616933169478728 x 10-5 —0.3616933215630078 x 107° 
2 y2 0.9815029948594308 0.9815030036954803 

YB 1.018493388207507 1.018493379371295 


i) 


0 0.5 1 15 


Figure 3: Variation of 2+ yi — ye — yg versus x which yi, yo and y3 are the 
numerical solutions obtained by the method of order 5 
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Figure 4: Variation of 2+ yi — ye — yg versus x which yi, yo and y3 are the 
numerical solutions obtained by the method of order 6 


The differential equations in Problem 82 satisfy a linear conservation law 
2+ yi(a) — yo(x) — s(x) = 0, (9) 


for all x. In Figure 3 and Figure 4, we have plotted the graph of 2+y1 —y2—y3 
versus x. We observe that for both methods of orders five and six equation (9) 
for the obtained numerical solutions holds approximately with high accuracy 
which demonstrate the accuracy of the applied methods. 


5 Conclusion 


For methods of higher orders (p > 5) with p = q¢ =r =, it is no longer pos- 
sible to solve the nonlinear systems of equations for satisfying RKS property 
by symbolic manipulation packages [5]. It seems that this difficulty does not 
appear for methods with fewer stages. In this paper we constructed RKS 
methods of orders p = 5 and p= 6 with r=s =3. 
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ee lo lE 5 sre Lb 
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awe gi 9 Cow bodes ale sl (SGLMs) e932 Gee eget bd cli-ss, Esl 2 ode 
BS- Sly Cole Cul chlo» or9 sluk — A oad cdl clea) cory cpl 09 ee 
dub gl jlade Sloe > ioly oad able cha, ohlS gol gh cl» wouc as ea ees 
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